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ABSTRACT 

We investigate the possible commensurabilities to be expected when two protoplanets 
in the Jovian mass range, gravitationally interacting with each other and an external 
protoplanetary disc, are driven by disc induced orbital migration of the outer proto- 
planct into a commensurability which is then maintained in subsequent evolution. We 
find that for a variety of protoplanet masses and typical protoplanetary disc proper- 
ties, as well as the setting up of 2:1 commensurabilities of the type recently observed 
in GJ876, 3:1 commensurabilities are often formed, in addition to 4:1, 5:1, and 5:2 
commensurabilities which occur less frequently. The higher order commensurabilites 
are favoured when either one of the planets is massive, or the inner planet begins with 
a significant orbital eccentricity. Detection of such commensurabilities could yield im- 
portant information relating to the operation of protoplanet disc interactions during 
and shortly post formation. 
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1 INTRODUCTION 

The recent discovery of a pair of extrasolar giant planets 
orbiting in a 2:1 commensurability around GJ876 (Marcy 
et al 2001) has raised interesting questions about the post- 
formation orbital evolution of this system. Assuming the 
system did not form in its currently observed state, the ex- 
istence of the commensurability indicates that disc induced 
orbital migration is likely to have occurred so as to gradu- 
ally reduce the planetary orbital separation until resonance 
was established. 

Simulations of two planets in the Jovian mass range in- 
teracting with a disc have been performed by Kley (2000) 
and Bryden et al (2000). The latter authors found a ten- 
dency for the two planets to open up gaps in their local vicin- 
ity, and for material between the two planets to be cleared, 
ending up interior to the inner planet orbit or exterior to 
the outer planet orbit, such that both planets orbit within 
an inner cavity. 

Snellgrove, Papaloizou & Nelson (2001) (hereafter pa- 
per I) performed a simulation of a system consisting of a 
primary star with two planets moving under their mutual 
gravitational attraction and forces produced by tidal inter- 
action with an externally orbiting gaseous disc. Angular mo- 
mentum exchange caused the outermost planet to migrate 
inwards until a 2:1 commensurability with the inner planet 
was reached. The subsequent dynamical interaction then re- 



sulted in the planets migrating inwards together maintaining 
the commensurability. 

The simulation results in paper I could be well matched 
with those of a simple A'^ body integration procedure used 
below, which incorporated simple prescriptions for the mi- 
gration and eccentricity damping of the outer planet due 
to interaction with the disc. Using this we investigate the 
range of commensurabilities that might be expected in two 
planet systems as a consequence of protoplanet disc inter- 
actions occurring in a standard protoplanetary disc model, 
where we assume a fixed migration time and consider various 
eccentricity damping rates. 

Depending on the nature of the interaction between the 
disc and the outer planet, as well as the details of the reso- 
nance into which the planets enter, we distinguish between 
four types of orbital evolution in the resonant migration 
phase which we denote as types A - D. We denote the or- 
bital elements of the outer planet with a subscript '1' and 
the inner planet with subscript '2'. Note that we find the 
evolution of the system depends on the details of the reso- 
nance into which the planets become locked. For higher or- 
der resonances such as 3:1 and 4:1, commensurabilities may 
be maintained in which differing resonant angles librate such 
that they do not cover the full (0, 2n) domain (see for ex- 
ample Murray & Dermott 1999 for an extensive discussion) . 
For a.p : q commensurability, we monitor the resonant angles 
defined by 
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Pp,q,k 



. = pXi — q\2 — pzoi + qzo2 + k{'uoi — ZU2 



(1) 



Here Ai,A2,ti7i and UJ2 denote the mean longtitudes and 
longtitudes of periapse for the planets 1 and 2 respectively. 
The positive integers p and q satisfy p > q, and there are 
p — q + 1 possible values of the positive integer k such that 
q < k < p. Although there are p — q + 1 corresponding an- 
gles, no more than two can be linearly independent. This 
means that if libration occurs, either all librate or only one 
librates. Both situations occur in our integrations. 
Type A: Eccentricities increase during migration until they 
reach steady state equilibrium values, with possibly small 
oscillations superposed. After this migration continues in a 
self-similar manner and all the resonant angles defined above 
are in libration. 

Type B: Eccentricities increase as the migration proceeds, 
until ei > 0.2 at which point we estimate the outer planet 
enters the outer disc. Our simple model breaks down at this 
point. Large values of ei arise when the eccentricity damp- 
ing rate is small. As in type A resonant migration all the 
resonant angles considered go into libration. The difference 
is that ei exceeds 0.2 before a steady state can be reached. 
An experimental model of the evolution of a two planet sys- 
tem corresponding to type B behaviour, and based on our 
knowledge of the migration of an eccentric planet interact- 
ing with a protoplanetary disc is presented in section 3.4. 
Type C: For higher initial values of 62, we find a mode of 
evolution in which ei and 62 rise continuously during the 
migration phase, even for efficient damping of eccentricity. 
The issue of the outer planet entering the disc again arises 
in type C migration. The primary difference between type 
C evolution and types A and B is that only one of the res- 
onant angles defined above is found to librate for the same 
apparent commensurability, so that the evolution differs. 
Type D: This mode of resonant migration corresponds to 
small values of ei < 0.2 being maintained for the outer 
planet, whilst the eccentricity of the inner planet grows to 
attain values close to unity. In this respect it differs from 
the other types of migration. This mode of evolution is re- 
lated to that described by Beust & Morbidelli (1996) when 
discussing the generation of star grazing comets in the (3 
Pictoris system. All the resonant angles librate in type D 
migration but with large amplitude. 



2 CONSERVATION OF ENERGY AND 
ANGULAR MOMENTUM 

We consider the consequences of a simple application of the 
conservation of energy and angular momentum. In the case 
when a near steady state for the orbital eccentricities is at- 
tained, a relationship between the orbital eccentricities and 
the circularization and migration rates induced by the disc 
is obtained for the case of type A migration. 

We consider two planets with masses m\,ni2, osculating 
semi-major axes ai,a2, and eccentricities 61,62. These orbit 
a central mass M*. The total angular momentum is 

J = J1 + J2 = mi^GM.ai{l - ei)+m2^GM,a2{l - el){2) 
and the energy E is 

GKUmi GM,m2 



We assume resonant self-similar migration in which 02/01, 
61, and 62 are constant. Then conservation of angular mo- 
mentum gives 



dJ _ J 1 dai I m2y/a2{l - e^) 
dt ^2ai dt \ miy/ai{l - e'f) 



-T 



and conservation of energy gives 

dE _ GALnii dai 
m ^ 2a| dt 



V mia2/ 



niT 



(4) 



(5) 



where ni = Ghh / a\ and we suppose there is a tidal 
torque —T produced by the disc which acts on mi. In ad- 
dition we suppose there to be an associated tidally induc ed 
orbital energy loss rate which is written as n\T / ^1 — e\ -\- 
D, with D = (GA/,mi6?)/(ai(l - el)tc). Here tc is the 
circularization time of mi that would apply if the tidal 
torque and energy loss rate acted on the orbit of mi with 
m2 being absent. In that case we would have dei/dt — 
—ei/tc- A migration time tmig can be defined through 
T — rriiy^ GM,ai{l — e1)/{3tmig)- This is the time for m 
to increase by a factor of e if m2 was absent and the eccen- 
tricity 61 was fixed. Note that tc and tmig are determined 
by the disc planet tidal interaction and may depend on 61. 
By eliminating ^ from (|) and (|) we can obtain a rela- 
tionship between 6i,62,tc, and tmig in the form 



el 



^1 -3/2 3/2 



i2 , fTllzfi)!^ 

^1 ^l-e2)l/2 



(6) 



We comment that this is general in that it depends only 
on the conservation laws and applies for any magnitude of 
eccentricity. However, it's derivation did assume self-similar 
migration with equilibrium eccentricities and so it does not 
apply in non equilibrium situations for which the eccentric- 
ities grow continuously. Note too that for a 2:1 commen- 
surability for which (ai/a2)'^''^ = 2, reduces in the case 
6^ <C 1 to the expression given in paper I which was obtained 
by analysis of the perturbation equations directly: 

2 tcm2ai 



3tmiff(2mia2 -I- m2ai) 



(7) 



E = 



2ai 



2a2 



(3) 



The above determines the eccentricity of the outer planet 61 
as a function of tc and tmio- 



3 NUMERICAL CALCULATIONS 

3.1 Model Assumptions and Physical Setup 

The basic assumptions of our model are that the two planets 
orbit within the inner cavity of a tidally truncated disc that 
lies exterior to the outer planet. Tidal interaction with this 
disc causes inwards migration of the outer planet on a time 
scale of tmig, and also leads to eccentricity damping of the 
outer planet on a time scale of tc . We assume that the inner 
edge of the disc is such that an eccentricity of ei > 0.2 will 
enable the outer planet to enter the disc, in basic agreement 
with the results of hydrodynamic simulations (e.g. Nelson et 
al 2000). 

We have performed three-body orbit integrations using 
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a fifth-order Runge-Kutta scheme, and the Burhsch-Stoer 
method as an independent check (e.g. Press et al. 1993). 
A torque was appUed to the outermost planet such that it 
migrated inwards on a time scale of tmig ~ 10* local or- 
bital periods and a damping force proportional to the ra- 
dial velocity was applied in the radial direction. This has 
a fixed constant of proportionality such that for small ei 
and in the absence of disc torques, the eccentricity damps 
on a time scale of tc — GOOA/" local orbital periods with val- 
ues of A/" = 1, 10, 100 having been considered. We comment 
that we have adopted the simplification of neglecting any 
possible dependence of the dam ping force and excepting the 
calculation presented in section 



3.4 



in 



on ei. 



This procedure was shown to be capable of matching 
the results of a detailed simulation in paper I, and the value 
of tmig adopted is consistent with that found from hydrody- 
namic simulations of protoplanetary discs interacting with 
protoplanets in the Jovian mass range (e.g. Nelson et al 
2000). Consideration of the orbital parameters of GJ876 (see 
paper I) suggests J\f ^ 1. However, there is some uncertainty 
in the functional dependence of the circularization rate of an 
orbiting protoplanet nonlinearly interacting with a tidally 
truncated protostellar disc as it would depend on the dis- 
tance to and form of the disc edge (Goldreich & Tremaine 
1980), and be sensitive to the presence of coorbital material 
which can lead to eccentricity damping through coorbital 
Lindblad torques (Artymowicz 1993). 

For length scale we adopt a fiducial radius, R, which 
for simplicity is taken to be 1 AU. However, scaling invari- 
ance allows this to be scaled to some other value if required. 
Runs were typically started with the outer planet at ra- 
dius 15.6-R while the inner planet was started at radius 
5R. A larger initial separation would be inconsistent with 
our model assumption of there being two planets orbiting 
within a tidally cleared cavity. Simulations by Bryden et 
al. (2000) suggest that the time scale to clear such a cav- 
ity becomes comparable to the migration time if the ratio 
of planetary semimajor axies becomes much larger than 3. 
The outer planet was assumed, because of interaction with 
the disc, to be in a circular orbit while starting eccentricities 
62 — 0.0, 0.05, 0.1, 0.2, 0.3 were adopted for the inner planet. 
For protoplanet masses we have considered all permutations 
of 0.4, 1 and 4 Jupiter masses assuming the central star mass 
is 1 Mq. 

3.2 Results 

The results of our survey are shown in table 1. In general 
for the shortest circularization rates with M — 1 and initial 
values of 62 < 0.1, self-similar (type A) migration occurred. 
For larger values of TV, type B migration was usually the 
preferred outcome. 

For larger initial values of 62 and — 1, either type 
A or type C migration occurred. For larger values of A/", 
types A, B, C, or D were all found to occur, with a strong 
preference for types B and C. Only one example of type D 
migration arose in our simulations. 
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Table 1. This table indicates the outcome and commensurabil- 
ity attained when a pair of masses mi , m2 measured in Jupiter 
masses becomes resonantly coupled due to disc driven migration 
of the outer planet. The initial eccentricity 62 and Af are indi- 
cated, ei is always initated as zero. The letters A, B, C, and D 
in columns 4-8 indicate whether types A, B, C, or D evolution 
occurred. 



Semlmojor Axis v's Time Eccentricity v's Time 
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Figure 1. This figure shows the evolution of the planet semima- 
jor axes and eccentricities for a pair of protoplanets with mi = 1, 
m2 = 1, initial 62 = 0, and A'' = 1. This case underwent self- 
similar (type A) migration in 2:1 resonance after having attained 
equilibrium eccentricities. The outer planet is denoted by the up- 
per line in the first panel and the lower line in the second panel. 



The unit of time is I0'^(i?/1AU) 



3/2 



yr- 



3.3 Initial 62 = 

For an inner planet on an initially circular orbit, a 2:1 com- 
mensurability normally resulted, except when one of the 
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I ^ 1 1 ^ 1 

500000 1e+06 1.5e+06 2e+06 2.5e+06 3e+06 

lime 

Figure 2. This figure siiows the evolution of the protoplanct 
semimajor axes and eccentricities for a pair of protoplanets, each 
of 1 Jupiter mass with J\f = 100 but with the migration rate 
of the outer planet modified such that it reverses for ei > 0.2. 
At small times, the upper two curves show logj^g(a/i?) while the 
lower two curves show ei, 62. The upper curve of the two denoting 
\o^lQ{a/ R) represents the outer planet, whereas at late times the 
lower curve of the two representing ei, 62 represents the outer 
planet. The unit of time is (_R/1AU)^/^ yr. 
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Figure 3. This figure shows the evolution of the protoplanct 
semimajor axes and eccentricities for a pair of protoplanets, with 
mi = 1, m2 = 4, initial 62 = 0.3, and M = 1. This case provides 
an example of type C migration in 3:1 resonance. Here the eccen- 
tricities of both planets grew continuously as the system evolved, 
even with significant damping of ei by the disc. Note that our 
simple model breaks down for e\ > 0.2. The outer planet is de- 
noted by the upper line in the first panel and the lower line in 
the second panel. The unit of time is 10"'' (R/IAU)^''^ yr. 

masses was 4 Jupiter masses in which case a 3:1 commensu- 
rability could occur. Plots of the evolution of the semimajor 
axes and eccentricities of a Jupiter mass pair with J\f = 1 
which undergo self-similar migration are illustrated in figure 
|l| In this case steady equilibrium eccentricities are obtained 
that are consistent with equations (^) and (^. We found 
that for the adopted migration rate, the transition between 
a 2:1 and a 3:1 commensurability occured for mi in the range 
2-3 Jupiter masses when 7712 was a Jupiter mass, with the 
transition mass being smaller for larger migration rates. 

3.4 Initial < 62 < 0.1 

When the inner planet was started with an initial eccentric- 
ity such that < e2 < 0.1, trapping in a 3:1 commensura- 
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Figure 4. This figure shows the evolution of the protoplanet 
semimajor axes and eccentricities for a pair of protoplanets, with 
mi = 4, m2 = 1, initial 62 = 0.3, and M = 10. This run pro- 
vides an example of type D migration in 4:1 resonance. Here the 
eccentricity of the inner planet grows without limit during the 
evolution, eventually reaching 62 ~ 1, whereas the eccentricity of 
the outer planet remains relatively small. As 62 approaches 1, a 
dynamical instability occurs, causing the planets to be scattered. 
The simplicity of our disc model does not allow this phase of 
evolution to be modeled accurately. The outer planet is denoted 
by the upper line in the first panel, and by the lower line in the 
second panel. The unit of time is 10^ (-R/IAU)'^''^ yr. 



bility was more common and could occur even for the lowest 
mass pairs. Also a few cases of trapping in a 4:1 commen- 
surability and one 5:2 commensurability were found. It was 
observed that the mode of migration in these cases usually 
corresponded to type A for A/" = 1, and type B for A/" = 10 
or 100. 

In order to investigate the potential outcome of a type B 
migration we considered a Jupiter mass pair with M = 100 
in which case the eccentricity of the outer planet grows to 
exceed 0.2. At this stage the assumed migration rate be- 
comes inapplicable. In fact migration may reverse on ac- 
count of penetration of the disc by the outer planet. In this 
context we note that a simulation of massive protoplanets 
by Papaloizou, Nelson & Masset (2001) indicated such a 
potential reversal for ei > 0.2 and linear torque calcula- 
tions by Paploizou & Larwood (2000) suggested migration 
reversal for an embedded protoplanet with eccentricity ex- 
ceeding a few times the disc aspect ratio. As an experiment 
we modified the disc induced migration rate of mi by a 
factor (1 — (ei/0.2)^). The resulting evolution is shown in 
figure H This adjustment can stall the migration rate as ei 
approaches 0.2 and lengthen the lifetime of the system be- 
fore a potential orbital instability occurs due to the devel- 
opment of large eccentricities in both planets. In fact in this 
case the eccentricity of the inner planet reached ~ 0.8 after 
~ 2 X I0''(ii/1AU)^/^ y. The system was found to be subse- 
quently stable for at least a similar time after a rapid disc 
removal. 

Thus we emphasise that the projected lifetimes of the 
resonantly migrating systems discussed here are dependent 
on the nature of the disc interaction, and at present are very 
uncertain. However, they could approach the protoplanetary 
disc lifetime if the mode of evolution discussed in the pre- 
ceding paragraph was to occur. 
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3.5 Initial 0.2 < 62 < 0.3 

For initial values of 62 in the range 0.2 < 62 < 0.3, it was 
observed that type C evolution became an important con- 
sideration, with only a single case of type D evolution being 
observed. For type C evolution, as in type B, the eccentrici- 
ties of both planets were observed to grow without reaching 
equilibrium, until ei > 0.2 at which point our simple model 
breaks down. The long term evolution in this situation may 
resemble that described in figure ^. We illustrate the be- 
haviour of type C evolution in figure ^. 

The single case of type D evolution obtained is illus- 
trated in figure ^ which should be compared with figure ^ 
since the only difference between these runs is the initial 
value of 62 (e2 = 0.2 in figure |^ and 62 = 0.3 in figure In 
this mode of evolution, the outer planet maintains a modest 
eccentricity while the eccentricity of the inner planet grows 
without limit. Although this mode of evolution appears to 
be rather rare based on the results of our calculations, it is 
interesting to note that it provides a method of generating 
very high eccentricities, such as that observed in the system 
HD80606 where the planetary eccentricity is e > 0.9 (Naef 
et al. 2001). 

Even for initially high values of 62, we still find a ten- 
dency for lower mass pairs of planets to enter lower order 
resonances such as 2:1, 3:1 and even 3:2, whilst the higher 
mass pairs tend to occupy the higher order commensura- 
bilities, even resulting in capture into 5:1 on a significant 
number of occassions. 



3.6 The Case of GJ876 

We have run a number of calculations to examine the pu- 
tative early evolution of the observed system GJ876, which 
exists in a 2:1 commensurability. The minimum masses of 
the two planets, in units where the central star s of solar 
mass, are mi = 6 Mjand m2 = 1.8 Mj. These parameters 
suggest that trapping in higher order resonance may have 
been expected rather than in 2:1. For initial conditions in 
which mi = 6, m2 = 1.8, initial 62 = 0, capture into 3:1 
resonance is favoured. For initial 62 ~ 0.1, 4:1 resonance re- 
sults. However, in a scenario in which the planetary masses 
were smaller during the initial resonant capture, capture into 
2:1 can occur. Simulations were performed with mi = 1, 
m2 — 1.8, initial 62 = and 0.1, which all favour capture 
into 2:1. Assuming a mass for the outer planet of mi — 3 
produces capture in 3:1. 

Two possible conclusions to draw from this are that: 

(1) . resonant capture occurred in this system when the plan- 
ets were of smaller mass, and that the mass of the outer 
planet increased due to accretion of gas from the disc. 

(2) . The planets formed sufficiently close together initially 
that capture into resonances of higher order than 2:1 was 
not possible. 



disc cavity while the disc interaction induces inward migra- 
tion. We find that as well as 2:1 commensurabilities being 
formed, 3:1 commensurabilities are obtained in addition to 
4:1, 5:1, and 5:2 commensurabilities which occur less fre- 
quently and usually for larger initial eccentricities of the 
inner planet. Possession of a small eccentricity ~ 0.1 before 
resonant locking takes place would favour 3:1 commensu- 
rabilities if migration over large separations has occurred. 
Initially near circular orbits would favour 2:1 commensura- 
bilities except when one of the planets is significantly more 
massive than one Jupiter mass. 

Up to now seven multiple planet systems have been de- 
tected with two showing a 2:1 commensurability (GJ876, 
IID82943 - see the website: 

obswww.unige.ch/~udry/planet/new_planet. html). It is too 
early at present to establish the frequency of occurence of 
commensurabilities in multiple systems but they may be 
fairly common. If orbital migration has occured over a large 
scale in extrasolar planetary systems the future detection 
of additional commensurabilities of the type discussed here 
is to be expected. We also comment that in some cases a 
high eccentricity for the inner planet orbit may be produced 
through resonant migration. If, depending on the details of 
disc dispersal and interaction, a scattering followed by ejec- 
tion of the outer planet occurs, a single planet remaining on 
an eccentric orbit may result. This will be a topic for future 
investigation. 
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4 SUMMARY AND DISCUSSION 

In this paper we have considered two protoplanets gravita- 
tionally interacting with each other and a protoplanetary 
disc. The two planets orbit interior to a tidally maintained 
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